Workflow for Ridgway’s Hawk Integrated Population Model
Brian W. Rolek
2024-07-31
Supplemental materials for: Rolek, B.W., McClure, CJW, Dunn, L., Curti, M., … Ridgway’s Hawk IPM and PVA
Contact information: rolek.brian@peregrinefund.org
Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/XXXXX.
The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/XXXXX/.
A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX
Plot model estimates
#load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_sites.rdata")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm1.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
out <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]),
as.mcmc(post[[5]]) ) #,
# as.mcmc(post[[6]]),
# as.mcmc(post[[7]]),
# as.mcmc(post[[8]]),
# as.mcmc(post[[9]]),
# as.mcmc(post[[10]]))
# Identify chains with NAs that
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
post2 <- post[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist## [1] TRUE TRUE TRUE TRUE TRUE# default settings for plots
plt <- function(object, params,...) {
MCMCplot(object=out,
params=params,
guide_axis=TRUE,
HPD=TRUE, ci=c(80, 95), horiz=FALSE,
#ylim=c(-10,10),
...)
}Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.
# Abundance of females at Los Haitises
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 1]"),
main="First-year (FY) minus hacked\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 1]"),
main="Adult nonbreeder (NB)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 1]"),
main="Adult breeder (B)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 1]"),
main="Adult Breeders and Nonbreeders\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 1]"),
main="All stages\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1]+datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
# Abundance of females at Punta Cana
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 2]"),
main="First-year (FY) plus hacked\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 2]"),
main="Adult nonbreeder (NB)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 2]"),
main="Adult breeder (B)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 2]"),
main="Adult Breeders and Nonbreeders\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 2]"),
main="All stages\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2]+datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
Population dynamics are determined by transitions, These transitions between stages are abbreviated as the starting life stage to the final life stage. For example a first-year recruiting to a breeder would be abbreviated as “FY to B”. I’ll list them here for convenience:
“FY to NB” is first-year to nonbreeder.
“NB to NB” is nonbreeder adult to nonbreeder adult.
“B to NB” is a breeding adult to a nonbreeder adult.
“FY to B” is first-year to breeder.
“NB to B” is nonbreeder adult to breeder adult.
“B to B” is breeder adult to breeder adult.
# Finer population segments
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 1]"),
main="Los Haitises\nFirst-years born",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY born",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
Other parameter estimates.
# I needed to abbreviate to save plot space
# FY=first-year, NB=Nonbreeder, B=Breeder
par(mfrow=c(1,2))
plt(object=out,
params=paste0("mus[",1:8, ", 1]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Los Haitises",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection")
)
plt(object=out,
params=paste0("mus[",1:8, ", 2]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Punta Cana",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
par(mfrow=c(1,1))
plt(object=out,
params="betas",
main= "Translocation effects",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
# Get predicted survival, recruitment, and detection
# with effects from tranlsocation and hacking
# Birds were only hacked from LHNP to PC here
# so we only predict values for PC
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in c(1,2,3,5,7)){
for (tr in 1:2){
pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] )
}}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
par(mar=c(8,5,1,1), mfrow=c(1,1))
for (tr in 1:2){
if (tr==1){
plot(1:5+c(-0.1,0.1)[tr], mus.md[c(1,2,3,5,7),tr],
pch=c(16, 17)[tr], cex=2,
ylim=c(0,1),
xlim=c(0.5, 5.5),
cex.axis=2, cex.lab=2,
ylab="Probability", xlab="",
main="",
xaxt="n", yaxt="n")}
else{
points(1:5+c(-0.1,0.1)[tr], mus.md[c(1,2,3,5,7),tr],
pch=c(16, 17)[tr], cex=2)
}}
axis(1, at=1:5, cex.axis=1.5,
labels=c("Survival\nFY", "Survival\nNB", "Survival\nB",
"Recruitment\nNB to B",
"Detection\nNB"), las=2)
axis(2, at=c(0, 0.25, 0.5, 0.75, 1), cex.axis=2,
labels=c(0, NA, 0.5, NA, 1))
for (m in 1:5){
for (tr in 1:2){
lines(rep(c(1:5)[m] + c(-0.1,0.1)[tr],2), mus.HDI95[1:2,c(1,2,3,5,7)[m],tr], lwd=3)
lines(rep(c(1:5)[m] + c(-0.1,0.1)[tr],2), mus.HDI80[1:2,c(1,2,3,5,7)[m],tr], lwd=6)
}}
legend(x=1.15,y=0.4, pch=c(16,17), pt.cex=2, cex=1.5, xpd=NA, horiz=T,
legend=c("Unmanaged", "Translocated" ))
# Fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("lmu.f"),
labels= c("Fecundity\n(log scale)\nLos Haitises",
"Fecundity\n(log scale)\nPunta Cana"))
f <- exp(outp$lmu.f)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)
# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.f))
for (s in 1:2){
f.pred[s,] <- exp(outp$lmu.f[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)
par(mar=c(6,5,1,1))
plot(c(1, 2)-0.1, f.md,
pch=16, cex=3,
ylim=c(min(f.HDI95), max(f2.HDI95)),
xlim=c(0.5, 2.5),
ylab="Fecundity", xlab="",
main="",
cex.axis=2, cex.lab=2,
xaxt="n", yaxt="n")
axis(1, at=c(1,2), cex.axis=2,
labels=c("Los Haitises\nNational Park","Punta Cana"),
padj=1)
axis(2, at=c(0, 0.2, 0.4, 0.6, 0.8), cex.axis=2,
labels=c(0, NA, 0.4, NA, 0.8))
lines(c(1,1)-0.1, f.HDI95[,1], lwd=3)
lines(c(2,2)-0.1, f.HDI95[,2], lwd=3)
lines(c(1,1)-0.1, f.HDI80[,1], lwd=6)
lines(c(2,2)-0.1, f.HDI80[,2], lwd=6)
points(c(1.1, 2.1), f2.md,
pch=17, cex=3)
lines(c(1,1) +0.1, f2.HDI95[,1], lwd=3)
lines(c(2,2) +0.1, f2.HDI95[,2], lwd=3)
lines(c(1,1) +0.1, f2.HDI80[,1], lwd=6)
lines(c(2,2) +0.1, f2.HDI80[,2], lwd=6)
legend(x=1.9,y=0.8, pch=c(16,17), pt.cex=3, cex=1.5,
legend=c("Untreated", "Treated" ) )
# Is fecundity at LHNP greater than PC
par(mfrow=c(1,1))
fdiff <- f[1,]-f[2,]
hist(fdiff, main="Fecundity difference\nbetween untreated sites")
# pd=
mean(fdiff>0)## [1] 0.9787# How big are the treatment effects on fecundity at each site
median(f.pred[1,]/f[1,])## [1] 6.147469median(f.pred[2,]/f[2,])## [1] 6.147469# gamma = nest treatment effect on fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("gamma"),
main="Anti-Parasitic Fly\nTreatment Effects", ylim=c(0,3))
par(mfrow=c(1,1))
sds <- paste0("sds[", 1:9, "]")
plt(object=out, params=sds,
exact=TRUE, ISB=FALSE,
main="Temporal SDs (synchrony among sites)",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection",
"Fecundity"))
sds2 <- paste0("sds2[", 1:9, "]")
plt(object=out, params=sds2,
exact=TRUE, ISB=FALSE,
main="Site-temporal SDs",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection",
"Fecundity"))
# Correlations among vital rates
# Plot is messy with only a few strong correlations
ind <- 1
Rs <- R2s <- c()
for (i in 1:(nrow(outp$R)-1)){
for (j in (i+1):nrow(outp$R)){
Rs[ind] <- paste0("R[",i,", ", j, "]")
R2s[ind] <- paste0("R2[",i,", ", j, "]")
ind <- ind+1
}}
par(mfrow=c(2,1))
plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time (synchrony)",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time (synchrony), continued...",
xlab = "Rhos", guide_lines=TRUE)
par(mfrow=c(2,1))
plt(object=out, params=R2s[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=R2s[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites, continued ...",
xlab = "Rhos", guide_lines=TRUE)
# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
main="First-year survival", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiA", ylim=c(0,1),
main="Adult nonbreeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiFYB", ylim=c(0,1),
main="First-year to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiAB", ylim=c(0,1),
main="Adult nonbreeder to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiBA", ylim=c(0,1),
main="Adult breeder to nonbreeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pA", ylim=c(0,1),
main="Nonbreeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.f",
main="", labels=labs,
xlab = "Year", ylab= "Fecundity")
abline(v=13.5, lwd=2)
Stage Structure of Each Population
mdFY <- apply(outp$NFY, c(1,2), median)
mdB <- apply(outp$NB, c(1,2), median)
mdF <- apply(outp$NF, c(1,2), median)
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number")
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
# Use median number of females in each stage
# to plot an approximate population structure
ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=Year)) +
geom_bar(position="fill", stat="identity") +
ylab("Proportion of population") +
facet_wrap("Site")
ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=Year)) +
geom_bar(position="stack", stat="identity") +
ylab("Numer of females") +
facet_wrap("Site", scales = "free_y")
Print parameter estimates
Parameter estimates for input into a population viability analysis.
pars1 <- c("sds", "sds2","mus", "betas",
"NFY", "NF", "NB", "Ntot",
"mn.phiFY","mn.phiA", "mn.phiB",
"mn.psiFYB", "mn.psiAB", "mn.psiBA",
"mn.pA", "mn.pB")
# Estimates for the survival model
# In this order: FY survival, NB survival, B survival,
# FY to B recruitment, NB to B recruitent, B to NB recruitment,
# Detection NB, Detection B
MCMCsummary(post2, params = c(sds, sds2),
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Mus are means for
# mus[1, site] , where site=1 is LH and site=2 is PC
# Survival of first years = mus[1,]
# Survival of nonbreeders = mus[2,]
# Survival of breeders = mus[3,]
# Breeding propensity of first years = mus[4,]
# Breeding propensity of nonbreeders = mus[5,]
# Transition from breeder to nonbreeder = mus[6,]
# Detection probability of nonbreeders = mus[7,]
# Detection probability of breeders
# modeled as logit(probability) = lmus[x,site] + betas[x]*translocated + eps[x,t] + eta[x,s,t]
# except breeder to nonbreeder recruitment = lmus[6,site]
MCMCsummary(post2, params = pars1[3],
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")MCMCsummary(post2, params = "lmus",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")MCMCsummary(post2, params = pars1[4],
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Fecundity
# modeled as log(f) = lmu.f[site] + gamma*treatment + eps[x,t] + eta[x,s,t]
# log scale
MCMCsummary(post2, params = "lmu.f",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Estimates of population size
# NFY= first year,
# NF = nonbreeders,
# NB=breeders,
# Ntot=total
# All are presented as N[time, site],
# where time 1=2011 ... 13=2023,
# site 1 = LH and site 2=PC
MCMCsummary(post2, params = pars1[5:8],
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")MCMCsummary(post2, params = pars1[9:16],
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Correlations among demographic rates time (synchrony)
MCMCsummary(post2, params = Rs,
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Correlations among demographic rates site x time
MCMCsummary(post2, params = R2s,
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")Model diagnostics
Check Goodness-of-fit
Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.
# Goodness of fit check
fit.check <- function(out, ratio=FALSE,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
jit=100,
ind=1,
lab=""){
par(mfrow=c(1,1))
# plot mean absolute percentage error
samps <- MCMCpstr(out, "all", type="chains")
rep <- samps[name.rep][[1]][ind,]
obs <- samps[name.obs][[1]][ind,]
mx <- max(c(rep, obs))
mn <- min(c(rep, obs))
plot(jitter(obs, amount=jit),
jitter(rep, amount=jit),
main=paste0("Mean absolute percentage error\n",lab),
ylab="Discrepancy replicate values",
xlab="Discrepancy observed values",
xlim=c(mn,mx), ylim=c(mn,mx),
pch=16, cex=0.5, col="gray10")
curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
bp1 <- round(mean(rep > obs),2)
loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
if (ratio==TRUE){
t.rep <- samps["tvm.rep"][[1]][ind,]
t.obs <- samps["tvm.obs"][[1]][ind,]
# plot variance/mean ratio
hist(t.rep, nclass=50,
xlab="variance/mean ", main=NA, axes=FALSE)
abline(v=t.obs, col="red")
axis(1); axis(2)
}
return(list('Bayesian p-value'=bp1))
}
# check goodness-of-fit for brood size
# breeder, ind=1
# fit.check(out, ratio=F,
# name.rep="dmape.rep",
# name.obs="dmape.obs",
# ind=1,
# lab="Breeder counts- Poisson", jit=300)
# # nonbreeder, ind=2
# fit.check(out, ratio=F,
# name.rep="dmape.rep",
# name.obs="dmape.obs",
# ind=2,
# lab="Nonbreeder counts- Poisson", jit=300)
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=1,
lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.04# first-year, ind=2
# poisson failed fit test bp=0
# Currently running models to try and fix
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=2,
lab="First-year counts\nNeg binomial-Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.38# fecundity
fit.check(out, ratio=F,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
ind=1,
lab="Fecundity-Neg binomial", jit=300)
## $`Bayesian p-value`
## [1] 0.81Examine Traceplots
Traceplots provide evidence that models adequately converged.
MCMCtrace(post2, pdf=FALSE, params= "sds")
MCMCtrace(post2, pdf=FALSE, params= "sds2")
MCMCtrace(post2, pdf=FALSE, params= "mus")
MCMCtrace(post2, pdf=FALSE, params= "betas")
MCMCtrace(post2, pdf=FALSE, params= "NF")
MCMCtrace(post2, pdf=FALSE, params= "NFY")
MCMCtrace(post2, pdf=FALSE, params= "NB")
MCMCtrace(post2, pdf=FALSE, params= "R")
